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Abstract 

Mono- and poly-disperse melts of oligomers (average length 10 monomers) of trans-1, 4-polyisoprene 
are simulated in full atomistic detail. The force-field is developed by means of a mixture of ab 
initio quantum-chemistry and an automatic generation of empirical parameters. Comparisons 
to NMR and scattering experiments validate the model. The local reorientation dynamics shows 
that for C— H vectors there is a two-stage process consisting of an initial decay and a late-stage 
decorrelation originating from overall reorientation. The atomistic model can be successfully 
mapped onto a simple model including only beads for the monomers with bond springs and 
bond angle potentials. End-bridging Monte Carlo as an equilibration stage and molecular dy- 
namics as the subsequent simulation method together prove to be a useful method for polymer 
simulations. 
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1 Introduction 



The understanding of polymer materials from the very local to the macroscopic scale is at the 
focus of theoretical material science. Simulations are a useful means for reaching this goal. Fully 
detailed simulations incorporating interaction centers for all atoms allow extensive investigations 
of polymer-specific models to compare directly to high resolution experiments like neutron scat- 
tering, 1 positronium annihilation spectroscopy 2 or nuclear magnetic resonance. 3 In contrast to 
simple generic models, 4 atomistic simulations capture the differences between polymer species. 
This understanding is a prerequisite for the design of tailormade materials for special applica- 
tions. 

The present contribution investigates, as an example, oligomers of trans-1, 4-polyisoprene. 
We focus on local structure and dynamics as the local scale is influenced mainly by the specific 
chemical structure. In this range all-atom models are necessary for a reliable static and dynamic 
investigation. 5 Moreover we compare our results to NMR measurements, which always sample 
atom-atom vectors. Therefore, it is useful for the analysis to have all atoms present. For the 
properties under investigation short chains are sufficient, as the connectivity and non-crossability 
contribute to phenomena on the generic large-scale level 6 ' 7 which is not of interest here. Several 
experiments have been performed aimed at the local properties of polyisoprene like reorienta- 
tion times and the packing. 8-13 Additionally, simulations of pure cis-polyisoprene are already 
known. 1,14 To our knowledge there are, however, not yet simulations of the trans conformer. 
Industrially this conformer is not as important as the ds-conformer. However, typical polyiso- 
prene materials include several percent trans content. In order to understand the influence of 
this content, it is worth investigating pure trans-p olyisoprene and look for differences from the 
other conformers. 

2 The simulation model and technical details 

2.1 Starting structures from quantum chemistry 

Three different chain length mixtures, each containing 100 oligomers (13200 atoms) of poly- 
isoprene with an average molecular weight of 682 g/mol (corresponding to 10 monomers) are 
simulated at temperatures of 300 K and 413 K. One system is mono-disperse (in the following 
referred to as system 1) and is simulated at constant density of 890 kg/m 3 (experimental density 
of the cis-conformer: 900 kg/m 3 ) 15 in an orthorhombic box under periodic boundary conditions 
with box lengths: 4.97nm x 5.16nm x 4.97nm. This is about two and a half times the end-to-end 
distance of the chains, thus artifacts from interactions of chains with their mirror images are 
absent. The torsional conformations were set up in the equilibrium distribution derived from 
quantum chemical calculations. Twenty-seven configurations of a dimer were investigated with 
quantum chemical calculations (for details see Section 2.3). Only seven of them are actually at 
low energy (Table 3) and only their statistics was used to generate initial torsional distributions. 
The chains were packed in a simulation box of 9nm x 9nm x 9nm and compressed to the final 
density in steps, using molecular dynamics. During this procedure anisotropic box fluctuations 
were allowed to better equilibrate stresses in the box leading to the orthorhombic box explained 
above. At the same time the time-step was increased from 0.01 fs to 1 fs. 

2.2 Starting structure via End-bridging Monte Carlo 

The two polydisperse systems (systems 2 and 3) were additionally equilibrated using the end- 
bridging Monte Carlo (EBMC) procedure 16 ' 17 before the molecular dynamics simulations were 
performed at a constant pressure of 101.3 kPa to determine the density in a poly-disperse melt. 
We used two systems with slightly different chain length distributions (see below) in order to 
increase the statistics and to look for influences of the actual realization of the ensemble. 

The EBMC was performed in a united-atom model as the positions of the hydrogens are not 
important at this stage. The non-bonded interaction potential parameters for this procedure are 



Interaction 



a nm 



e[kJ/mol] 



C-C,C-CH,CH-CH 
C-CH 2 ,CH-CH 2 
C — CH3, CH — CH3 



0.38 

0.4257 

0.4257 

0.45 

0.45 

0.45 



0.4186 
0.4249 
0.6299 
0.3918 
0.6095 
0.9479 



CH2— CH2 
CH2 — CH3 
CH 3 — CH 3 



Table 1: Nonbonded parameters (Lennard- Jones 12-6) used in the united atom end-bridging Monte 
Carlo simulations of irans-polyisoprene. Bond lengths, angles and torsion potentials are the same 
as in the all atom simulations. All interactions within a topological distance of up to 3 bonds are 
excluded, as they are taken into account in the torsion potentials. 

shown in Table 1. All Monte Carlo steps are accepted according to the Metropolis criterion. 

The end-bridging steps were performed using a hypothetical chain (see dashed line in Fig- 
ure 2). The monomers on this chain consist only of the carbons 1, 4, and 5. Unlike the original 
procedure for polyethylene, 16 ' 17 moves can only be performed breaking the bonds between car- 
bon 5 and 1 as the topology of trans-p olyisoprene must not be altered. As the bridging procedure 
always needs three trimers: one to bridge from, one to bridge to, and one as the bridge, the three 
carbons of one pseudo-monomer act as such trimers. After the respective move the positions of 
carbons 2 and 3 have to be recalculated. Since the double bond keeps the whole monomer in 
plane the positions are defined exactly if bond lengths and angles are left unchanged. In the 
Monte Carlo procedure no bond length or bond angle was changed, neither for the hypotheti- 
cal chain, nor for the atomistic chain. Any of these local degrees of freedom were left for the 
molecular dynamics to equilibrate. 

Anisotropic box fluctuations were allowed but no shearing. The systems are slightly poly- 
disperse, as this version of EBMC does not work at constant topology, i.e. chain length distri- 
bution. The molecular weight distribution is sharply peaked at A = 10 (Figure 1) as we did 
not intend to equilibrate the chain length distribution but to investigate nearly mono-disperse 
systems. We started with a monodisperse sample and let the end-bridging procedure run for 
a limited time with chemical potential fi = for 5 < N < 15 and \x = — 00 elsewhere. Strict 
monodispersity, however, would require much longer equilibration times as EBMC would be 
prohibited. On average, every chain underwent 2.8 successful end-bridging moves in addition to 
many local moves before the molecular dynamics started. The resulting distribution shows two 
weak side-peaks at N ss 6 and N ~ 14. These can be explained by the procedure. Chain lengths 
of 9 or 11 are not directly accessible by endbridging moves starting from the initial monodis- 
perse (N = 10) sample, as no single monomers can be transferred. Since they have to be reached 
indirectly, they have a low probability. 



As the torsion potentials are very important for the configurations of the chains, quantum chemi- 
cal calculations were performed with the packages Gaussian 94 18 and Gaussian 98. 19 The energy 
differences between different minima were calculated with hybrid density functional calculations 
(B3LYP) 20 and the barrier heights are calculated by Hartree-Fock using a 6-31G** basis set on 
a dimer of trans-polyisoprene. In order to calculate barrier heights, constrained optimizations at 
fixed bond lengths are applied where the torsions are changed in steps of 15 or 30 degrees. The 
results were fitted to a Fourier series in the torsion angle with four terms. The constant energy 



2.3 The force-field 
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Figure 1: Molecular weight distributions of systems 2 (filled bars) and 3 (open bars) 
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Table 2: Force-field parameters for the torsion angles. For the nomenclature see Figure 2. 



shift is discarded as it does not enter the forces 



Vtors = ^2 k T /2 1 - cos[i(r - T )] 



i=l 



(1) 



The resulting parameters are shown in Table 2. The relevant minima are first supposed to occur 
in trans and gauche states. The initial torsion states are set up according to this distribution 
(before EBMC) in a Markov chain way, i.e. the chains are built monomer by monomer and the 
triple of torsions connecting two monomers is selected by random numbers according the energy 
distribution. 

The bond angle potential and the improper dihedral potential, which keeps the atoms at the 
double bond in plane, are harmonic with force constants adapted from previous simulations of 
small molecules (Table 4). 21,22 The harmonic dihedral has a strength of 160 kJ/(mol*rad 2 ) and is 
applied to the dihedrals Ci— C2— C3— C4, Hc2~ C2— C3— C5, and C2— C3— C4— C5. Also the non- 
bonded interaction parameters (Table 5) are adapted from simulations of cyclohexene. 22 They 
were derived using the automatic simplex parameterization. 21 The polymer density was not fully 
satisfactory with the original cyclohexene force-field. 23 Thus, the a value of the hydrogens was 
slightly increased to reproduce the experimental density at 300 K. Lennard-Jones interactions 
between unlike atoms were based on the Lorentz-Berthelot mixing rules. 24 

Atoms connected by any bonding potential did not interact by the Lennard-Jones potential. 
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1 
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180.0 
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0.33 
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133.5 
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0.00 
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-70.5 
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4.91 
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-86.0 


-59.0 


-108.3 


6.72 



Table 3: Relevant torsion conformations after geometry optimization using the hybrid method 
6-311G**/B3LYP. 




Figure 2: Carbons in the dimer of irans-polyisoprene with definition of the torsions. Additionally, 
the hypothetical chain for the end-bridging moves is shown as dashed line. 
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Table 4: Angles and bond lengths for atomistic simulations. The equilibrium values of the angles 
0o and the bond lengths l b are derived from density functional calculations using Gaussian 94 18 and 
experimental data, 46 k is the force constant for the harmonic bond angle potential. The bond lengths 
are constrained. 
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Table 5: Force-field parameters for the non-bonded interactions, m is the atom mass, a the interaction 
radius, and e the interaction strength. 
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300 


1184 


1.00 


890 
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300 


2012 


1.05 


917.4 
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300 


1737 


1.05 


916.8 


3 


413 


792 


1.05 


826 



Table 6: Simulation parameters for the three systems under study. Before the above-mentioned 
simulation time was started a few hundred picoseconds equilibration time was waited for. 

Additionally, the following non-bonded interactions were excluded: all within one monomer, and 
all C— C, C— H and H— H interactions between the second half of the carbons of one monomer 
(atoms C3, C4, and C5) and the first half of its following neighbor (Ci and C2) including the 
hydrogens connected to them for system 1. For system 2 and 3 the latter C— H and H— H 
interactions were not excluded (only up to "1 — 4" interactions), which leads to differences as 
the different torsion states alter the distances between these atoms so that the effective energies 
of the torsions are shifted. This happened due to an oversight in the initial simulations. As the 
simulations are very time consuming we deemed a complete repetition not necessary. 

Constant temperature and pressure are ensured using Berendsen's method 25 with time con- 
stants 0.2 ps for temperature and 8 ps (system 2) or 12 ps (system 3) for pressure, respectively. 
The pressure coupling using a compressibility of 2 x 10~ 7 kPa -1 was employed for the three 
directions independently. All simulations were performed using the YASP molecular simula- 
tion package 26 with a time-step of 1 fs and a cutoff for the non-bonded interactions at 0.9 nm. 
Configurations were saved every picosecond. Simulation lengths and polydispersities are shown 
in Table 6. Bond lengths are constrained to the desired values of Table 4 using the SHAKE 
algorithm. 27,28 

3 Thermodynamic and static structural properties 
3.1 Density 

The systems simulated under constant pressure conditions arrived at densities of 917.4 kg/mol 
(system 2) and 916.8 kg/mol (system 3), respectively, which is a discrepancy of less than 2% 
relative to the experimental value of 900 kg/mol for chains of 16 monomer length. 15 The isochoric 
simulation of system 1 had a pressure of -500 kPa. In NVT simulations, a negative pressure of 
this magnitude means that the system is not exactly at the correct density but would like to 
contract a little further. As all densities are quite close to the experimental value, which itself is 
not too certain, since it was determined for a mixture of cis- and trans-PI, a closer refinement 
of the force-field parameters was not deemed necessary. The density also depends weakly on the 
intra-chain part of the potential. If the torsion potentials are switched off, the density increases 
by about 2%, as the chains can pack more effectively (System 1). 




Figure 3: Bond vector correlation functions in the atomistic polyisoprene simulations of system 1. The 
filled symbols are correlations between direction vectors in the monomers, (diamond: double-bonds, 
circle Ci— C5) The open symbols correspond to vectors connecting atoms belonging to neighboring 
monomers (squares: C2, triangles Ci). 

3.2 Single chain properties 

The mean squared end-to-end distance of the oligomers was measured between the terminating 
carbons (C™ 0110 1 and C™ ono n ). The respective result of 5.12 nm 2 corresponds to 0.0075 nm 2 mol/g 
for the monomer- weight-specific end-to-end distance. The experimental value of 0.0060 nm 2 mol/g 
per monomer is for a mixture of trans and cis polyisoprene. The experiments are performed doing 
small angle neutron scattering in a melt of 7-mers. 29 

The persistence length l p measures direction correlations of unit vectors along the chain. It 
is calculated using suitably defined points along the backbone. The function 

(u(r) ■ u(fo)) = e 1 p (2) 

is fitted against the bond correlation. There is some freedom in how to define the tangent 
vector in an atomistic model. The vectors connecting the Ci (or the C2) atoms of adjacent 
monomers were investigated. Additionally, intra-monomer vectors along the double-bond and 
vectors connecting the two end carbons of the same monomer (Ci— C5) are included. The bond- 
correlation functions are not really exponential (Figure 3), as the interactions are complex and 
the chains are short. Thus, persistence lengths deduced from the fitting procedure can only 
be taken as estimates. End-effects were not excluded for reasons of statistics. The correlation 
functions are different for different vectors. The vector representing the shortest connection (the 
double-bond) has the shortest correlation length. The persistence lengths range between 0.5 (for 
the double bonds) and 1.5 (for the intermonomer vectors) in monomer diameters (~ 0.3 nm). 
Polyisoprene is, therefore, rather flexible already on the length scale of a few monomers. The 
monomer itself is intrinsically stiff, but the three torsions between neighboring monomers provide 
a flexible link. 



3.3 Local packing of chains 

The local structure in the melt can be characterized by different pair distribution functions. 
Figure 4 shows inter-chain radial distribution functions of the different atom pairs present in 
polyisoprene. The curves between the different systems differ only slightly due to density differ- 
ences. The conspicuous absence of a distinctive first peak in the all-atom and the hydrogen g(r) 




Figure 4: Interchain radial distributions of the different atomic pairs (full line system 1, broken line 
system 2) a) all atoms, b) only carbons, c) only carbon hydrogen d) only hydrogen pairs 
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Figure 5: Center of mass radial distribution function of the atomistic polyisoprene chains. A corre- 
lation hole on local scales is seen, on bigger length scales the distribution is flat. The dashed line 
corresponds to system 1, the solid line to system 2, the dotted line to system 3. 



at short distances indicates that chains cannot easily interpenetrate. Structure is averaged out 
in the all-atom g(r), in contrast to, e.g., the carbon-carbon distribution function, which exhibits 
several clear peaks. However, also these do not rise to values gcc( r ) > 1- 

Figure 5 shows the center-of-mass distribution of whole chains. In spite of the statistical noise 
it is clear that they can approach as close as 0.2 nm. The similarity of the RDFs for systems 2 
and 3 suggests that the combination of EBMC and MD leads to a reliable structure. At distances 
of more than 1 nm the distribution is quite flat, at shorter distances a correlation hole is clearly 
visible. System 1, which was not equilibrated using the end-bridging, exhibits unrealistically 
sharp peaks. The differences in the molecular weight distribution also contribute to the different 
structures. However, molecular dynamics alone is not able to equilibrate a simulation of this 
size. The introduction of EBMC brings us a good step further, although some remnants of the 
setup may still be present. 

Carbon— carbon radial distribution functions (RDF) resolved according to carbon type al- 
low the study of preferential arrangements between different chains. Partial pair distribution 
functions between the five different carbons present in polyisoprene were recorded (Figure 6). 
The methyl carbon (C4) is the most exposed and can, therefore, approach closest to the others 

0.4nm). At this distance there is also a shoulder in the Ci, C2 and C5-RDF indicating direct 
contact. C3 is the most "shielded" carbon with a slight shoulder at direct contact to Ci. It 
is linked to C2, C4 and C5, thus, it is often found as second contact (sa 0.55nm). The two 
methylene carbons Ci and C5 are coordinated very similarly. C2 is easily accessible, as it has 
only one hydrogen but not very exposed to other chains, since it is part of the double bond in 
the backbone leading to weak structure. 

Integration of the pair distribution function yields the number of neighbors of an atom in 
a shell of radius r. By relating the local number density pi oc &\ = 4/3^,3 1 n c being the number 
of carbons, to the overall concentration p(oo) = 8.102nm -3 , local enrichment (x := ^'gyj > 1) 
and depletion (x < 1) can be resolved. Overall integrated values in Figure 6f are smaller than 
unity due to the correlation hole, as only foreign chains are taken into account. In the innermost 
shell (r < 0.45nm, cf. Figure 6f) of all carbons methyl groups (C4) are enriched, sp 2 carbons 
(C2 and esp. C3) are depleted, whereas the methylenes (Ci and C5) are close in concentration. 




Figure 6: Partial inter-chain radial distribution functions (system 2, T=300 K). a) Ci b) C2 c) C3 
d) C 4 e) C 5 . In all figures the definition of line styles in figure c) applies. Figure f (inset table): 
Integrated number of neighbors in the innermost shell r < 0.45nm, only foreign chains. 



In the second shell (0.45 < r < 0.65 nm), this is partly reversed, as C3 is enriched and C4 is 
weakly depleted. At distances larger than 0.65 nm all species are found at bulk concentration. 
In brief, one can say that the methylene carbons occur at constant density almost everywhere. 
The methyl and the double bonded carbons (esp. C3) show much more structure. Monomers 
of different chains thus approach each other typically with their side groups as closest contact. 
Orientational influences from the double-bond keeping the monomer planar play a role, too (see 
below). The whole overall structure in the RDFs extends about two monomer sizes (r < 1 nm), 
whereas the concentrations of different carbons level out already at 0.65 nm. 

The local structure is not fully described by the (spherically averaged) radial distribution 
functions. Mutual orientation is measured by the orientation correlation function OCF (Figure 7) 



The distance r is measured between the centers of mass of the respective pairs. Again there 
are several possible choices for the tangent vectors Uj. Orientation correlations of the vector 
connecting the methylenes (Ci— C5) extend over several inter-atomic distances (Figure 7a). They 
resemble the packing of model chains consisting of simple beads. 30 The very few segments that 
approach closely (cf. Figure 4 and 6) show a perpendicular orientation. A parallel ordering peak 
is encountered at about 0.4 nm. The first Legendre polynomial (-Pi(r) := (tti ■ Uj), Figure 7b) 
which carries direction information shows that these inter-chain contacts have a small preference 
of neighboring chains running in the same direction. At larger distances the explicit atomistic 
structure is no longer important and the structures become much broader. But there are still 
structural effects originating from the packing visible up to about 0.7 nm, about two chain 
diameters. 

The order of the double bonds between the chains exhibits more atomistic details. As the 
double bond lead to a planar configuration of the environment, there is parallel orientation of 
the neighbors especially at short distances (r < 0.4 nm) which is visible in the first Legendre 
polynomial. The bond vectors C5— Ci (between subsequent monomers), which are about the 
same length as the double bonds, show even less structure than the simple model. 

The inter- monomer vectors (Figure 7c), in contrast, display even more generic features which 
can also be seen in a simple model 3 ' 30 ' 31 (compare Section 5). As they describe larger segments 
compared to the vectors discussed above, the features are less accentuated. On close contact 
(r < 0.4nm) an almost perfect perpendicular order is found. At intermediate distances (0.4nm < 
r < 0.5nm) a preferred parallel alignment, and for the case of the C2 — C2 vector a second 
perpendicular region appears (r ~ 0.8nm). The differences between the two different inter- 
monomer (Ci— Ci and C2— C2) vectors are small. So for orientations on length scales as small 
as monomer sizes, simple models already provide a useful description if the persistence lengths 
are the same. 

This shows that the generic packing effects are important for the understanding of the struc- 
ture of atomistic models. However, the fine structure at short distances, as found in the first 
Legendre polynomial, cannot be explained by generic arguments as here the detailed chemistry 
of the polymer is important. 

3.4 Structure of the melt 

Radial distribution functions are quite illustrative in characterizing the structure. However, the 
experimental observable in scattering experiments is the structure function. The static melt 
structure factor 



is shown in Figure 8 in comparison to simulations of cis-polyisoprene and experiments on a 
mixture dominated by the cis-conformer. 10 The melt structure factor shows a clear peak at 
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Figure 7: Static inter chain orientation correlations OCF (T=300 K, system 1). a) Double-bonds and 
vectors connecting C 5 to Ci of the next monomer between neighboring chains. Dot-dashed: atomistic 
double bond OCF; solid: C5— Ci vectors; dashed: OCF of a simple model with persistence length 
l p — 1.5 monomer units 31 scaled for coincidence at the first maximum to the solid one. b) Same as 
figure a) but Pi in order to show the directionality of the correlations, c) Inter-monomer vectors, solid 
line: Ci— Ci, dot-dashed line: C 2 — C 2 , dashed line: simple model for comparison. Note the different 
ordinate scale. 
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Figure 8: a) Static structure factor of trans polyisoprene from this work (solid line: T=300 K, dashed 
line T=413 K) with the scattering lengths of all atoms taken to be the same in comparison to cis-PI 
(dotted line, T=413 K, data from ref. 1). b) Experimental structure factor obtained by neutron 
scattering at an unknown temperature at molecular weight of 17.200 (data from ref. 10). 



about 15 nm . In addition, there is some smaller structure, especially higher order peaks. The 
lower limit of resolution is given by the size of the box corresponding to a minimum | A; | -vector of 
about 0.47T nm -1 . At higher temperature the overall structure flattens out with less pronounced 
peaks, as expected. 

Moe and Ediger performed simulations on pure cis-polyisoprene at 363 K and 413 K with one 
chain of 100 monomers. 1 ' 14 This, obviously, reduces the influence of end effects. The densities 
were much lower than the experimental values (798 vs. 869 and 775 vs. 836 kg/m 3 , respectively). 
The structure functions of both simulations are, however, comparable at 413 K. Neither simu- 
lation compares really satisfactorily to the experimental structure factor. The double-maximum 
structure of the first peak is not reproduced, the lower maximum (~13 nm' 1 ) being enhanced, 
the higher (sa 18 nm -1 ) being reduced to a shoulder. The positions of the peaks, however, are 
in reasonable agreement. As the experimental temperature is not given, one cannot say whether 
at least part of the discrepancy is a low temperature effect, indicating the formation of a glass, 
whereas both simulations are performed deep in the melt. Another reason for the disagreement 
could be the fact that the short chains presented are still oligomers. The alternative would be 
one long chain, as in refs. 14 and 1. However, the periodic boundary replication of a single chain 
this imposes a strong artificial periodicity to a amorphous melt, which makes it difficult to look 
at inter-chain effects. 



4 Dynamical properties 



4.1 End-to-end vector 



Figure 9 shows for the two polydisperse systems 2 and 3 the reorientation correlation functions 
(first and second Legendre polynomial) 




(5) 



(6) 



of the end-to-end vector of the chains of length N = 10, which is defined as the vector connecting 
the two terminal carbons C™ 0110 1 and C" lono n . The relaxation time is clearly longer than the time 
accessible in the simulations. Even system 2, which was simulated for more than 2 ns, did not 




Figure 9: Relaxation of the end-to-end vector. Solid curves are first Legendre polynomials, dashed 
curves second Legendre polynomials. 

a) T=300 K: Of both pairs, the upper curve is the correlation function for system 3 and the lower 
one for system 2. 

b) For T=413 K (system 3) the relaxation times are much shorter. 



relax appreciably on this time scale; local vectors relax, of course, much faster. At the higher 
temperature of 413 K the relaxation time decreases drastically. Still, one has to be cautious 
discussing length scales of more than a monomer. This figure shows some scatter between the 
two systems, which may be taken as a rough estimate of the error of the simulations. 



4.2 C— H bond reorientation 

Reorientation in polyisoprene melts with 92% cis-conformer was measured by the group of 
Laupretre and Monnerie. 8,11 Another investigation with a higher trans content of 22% 9 focused 
also on the ds-conformer. Experimentally, the direct observable is the T\ time. 

i = ^tct| r J(wH _ wc) + 3J(u;c) + 6J( ^ h + wc) i _ (7) 

Ti Wr 2 c _ R L J 

The 7j are the gyro-magnetic ratios of the respective nuclei and cue and wh are the Larmor 
frequencies, r is the distance between the nuclei. The function J(u>) is the spectral density, i.e. 
the Fourier transform of C rc0 r of the respective C— H vector 

1 r°° 

J(oj) = - / C rcor (t)e wt dt . (8) 

In atomistic simulations, T\ has been determined for different polymers. 14,32 ' 33 Moe and 
Ediger use the limit of extreme narrowing (ujT rcor < 1) to analyze their cis-polyisoprene data 
at T = 413 K, as is done in most other atomistic simulations. 14 This has the advantage that 
T\ becomes independent of the Larmor frequencies, 34 which cannot be measured in simulations 
without extrapolation. For long chains in a simple model one sees that extrapolation starting at 
such high values is very questionable. 31 The spectral density J(lo) is for very short reorientation 
times independent of uj: 

± ' 1 rmr 



However, the extreme narrowing regime is not reached normally by the experiments, as high 
temperatures are needed to yield correlation times that the results correspond to the limit. 




Figure 10: Reorientation of C— H vectors in inner monomers of atomistic polyisoprene chains (sys- 
tem 3, chains of length 10,) a) Vinyl C2— H depending on monomer position (1: end monomer, 2: 
next to end monomer, etc.), b, c) Methylene groups of the central monomer at different temperatures. 



The reorientation time T re0 r is defined by the time integral over the correlation function 

/•oo 

7~reor — / Creordi , (10) 
JO 

and in the extreme narrowing limit this is directly linked to the T\ time for a C— H vector 

Tf 1 = 10nK Tlcor , (11) 

where K is a constant related to the bond length and n is the number of protons connected to 
the respective 13 C. Note that a shorter T\ corresponds to a longer reorientation time. 

The C— H vector reorientation is followed in the simulations. To minimize chain-end effects 
only the inner-chain monomers are included in Figure 10 b and c. Comparing Figures 12 a and 
10 a, one sees that the hydrogen connected to the backbone at carbon C2 is strongly tied to 
the backbone. Thus, even such local quantities as bond vectors can be used as observables to 
monitor the dynamics of intermediate-size chain segments. 

Figure 10 b illustrates that different torsion potentials and the side group have considerable 
influence on the reorientation of vectors. The reorientations of the two methylene carbons Ci 
and C5 differ. As C5 is vicinal to the methyl group, the steric hindrance for the hydrogens tied 
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Table 7: Comparison of the experimental (ds-PI) and simulation (cis and trans-PI, system 2) data 
for the efficiency of the initial stage of the reorientation process. a^ t originates from an exponential 
fit of the second stage extrapolated to t — and ai ps is the value of 1 — C rcor at lps. In the analysis 
of the simulations for ds-polyisoprene a stretched exponential second process was assumed. The 
experiments used a range of temperature between 283 K and 363 K (ref. 8). The trans simulations 
were at 300 K (this work) and the cis simulations at 363 K (ref. 1). 

to it is more pronounced. Thus, the initial stage of the reorientation, present for the Ci— H 
vector, is considerably smaller for C5— H. This is still visible at 413 K (Figure 10 c), although 
all relaxations are faster. The above-mentioned experiments 8 ' 9 on cis-polyisoprene melts show 
the same tendenies. 

Similar to our reorientation correlation function (cf. Figure 10) the experimental data 8 
provide evidence for a two stage process, the first part being simply exponential followed by a 
non-exponential long-term stage, described by model correlation function 

C rcor (t) = ae- tT0 + (1 - a)e- i / T2 e-*/ T1 / (t/r 1 ) , (12) 

with To the local libration time, T\ the time of conformation jumps, and Ti connected to damping; 
Iq is a Bessel function. The separation of time scales for polyisoprene is t±/tq > 150 for the two 
faster characteristic times; 11 the two slow processes (t±, T2) differ by a factor of 40. A separation 
of motions was used by Lipari and Szabo as well to analyze NMR data of polymers. 35,36 The 
correlation function they used has the simpler double exponential shape 

C Icor (t) = S 2 e~ t '^ + (1 - S 2 )e- 1 ^ , (13) 

where the generalized order parameter S is related to the parameter a of Laupretre et. al. 11 
The reorientation motion is described by one local and one global reorientation if there is no 
reptation. 

Actual numbers for the three times of the model (Eq. 12) are not provided in refs. 8, 11, 
but only values for a. Although the experimental polymer is mainly cis the values for a for the 
different vectors are comparable between our simulations and experiment (Table 7). The simu- 
lation data in our case were determined using a fit to Eq. 13 with a purely exponential form of 
the second term disregarding the Bessel function (Iq = 1). The second column in Table 7 shows 
the values of a estimated from the value of C rcor at 1 ps, which is the shortest time resolved in 
the simulations (ai ps = 1 — C rcor (lps)), coming closer to experimental values. This assumes that 
the first process is too fast to be resolved here. If the conformational jump time is regarded as 
the time for torsion rearrangement (see below), and keeping in mind that the two times differ 
at least by a factor of 150, the guess for ai ps is probably more realistic. Still, the simulations 
underestimate the difference between Ci and C2 and overestimate the one between Ci and C5. 
The data provided here correspond to a sample of pure trans-p olyisoprene oligomers which is pre- 
sumably quite different from real cis-p olyisoprene with some added irans-conformer. Moreover, 
the discrepancy becomes weaker for system 1 (below). Recent simulations on cis-polyisoprene 
at higher temperature (T=363 K and T=413 K) were interpreted in terms of a two stage corre- 
lation, too. There a separation between the exponential first stage and a stretched exponential 
second process was deduced. 1 The corresponding a-parameters are included in Table 7. Except 
for the C5— H vector they are comparable to our data. For C5— H they are even farther from the 
experimental value. 
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Figure 11: a) Short time reorientation of methyl C— H vectors, monomer number as in legend (end 
monomer vs. central monomer, System 3, T=300 K, only 10-mers are included), b) Long time. 
Clearly, the two stages are separate. A slight difference between the two terminal monomers is visible 



In the reorientation of the C4— H vectors of the methyl side group the two-stage reorienta- 
tion is clearly visible (Figure 11 a). The first process is non-exponential on the scale of a few 
pico-seconds. The monomer index has almost no influence, only the very local surroundings 
contribute. On this time scale, the vector does not experience the connection to the chain. 

The long time tail, however, is linked to the overall reorientation. The bonding to the chain 
leads to a bias in the orientation of the methyl group, which prevents total decorrelation on the 
short time scale. The second (exponential) process is influenced by chain end effects with decay 
constants of r cn d ~ 1 ns and r cen t C r ~ 3 ns, respectively (Figure 11 b). The decay times were 
determined by an exponential fit of the time region between 500 and 1250 ps. 

Correlation times for the reorientation of various C— H bonds in cis-polyisoprene melts were 
calculated by Moe and Ediger, who used one long chain under periodic boundary conditions. 14 
Table 8 compares the results of this work to their data and the extreme narrowing limit of 
the experiments. At 300 K, the correlation times were determined by numerical integration of 
C reor (t) over the first nano-second and an analytical correction for the exponential tail. For the 
methyl groups the numerical integration extended only to 20 ps. 

Heating to 413 K speeds up the simulation. Here the numerical integration was performed 
up to 100 ps, again with an analytical correction, except for the methyl groups (20 ps) and the 
Ci— H (200 ps) (label A in Table 8). For the oligomers, however, the results are very similar to a 
numerical integration to 400 ps, which was performed additionally in order to compare directly 
to the data by Moe and Ediger (label B in Table 8). 

The discrepancies between the different systems give a measure of uncertainty in the results. 
The reader is reminded that the systems are governed by slightly differing force-fields and are 
equilibrated in different manners. 

If the data is compared directly to extrapolated experimental data an overall discrepancy 
of about 50% is found between the trans-simulations presented here and the experiments. The 
integration error in the simulation as well as the extrapolation of the experiments are sources of 
error. The systems are not the same and the simulation model does not reflect reality perfectly. 
The experiments themselves are not perfectly reliable. Witt et al. showed that NMR experiments 
for systems as simple as liquid benzene can result in reorientation times differing by an order of 
magnitude. 37 

In the simulation system 2, the hydrogen at C5 reflects too much how the backbone reorients 
on the time scale of the double bond. Experimentally there is a difference between C5-H and 
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Table 8: Comparison of simulated reorientation times in ^rans-polyisoprene (this work) and cis- 
polyisoprene (ref. 14) to data determined by extrapolation of experiments at lower temperatures into 
the extreme narrowing limit. The cis- values 14 were determined by numerical integration of the first 
400 ps. For the £rans-polyisoprene only the two innermost monomers are used. At 413 K the analysis 
was done in two ways in order to compare more directly to the Moe and Ediger simulations. A: 
numerical integration and exponential long-time tail (see text), B: numerical integration to 400 ps. 
The simulation errors are estimated to be about 20% (difference between systems). The reorientation 
of the methyl group in the NVT simulation could not be calculated meaningfully, as there was a 
force-field problem. (The hydrogens were connected to the carbon with an additional torsion, which 
was too strong.) 

C2-H. This may result from the interaction with the methyl group which repels the H at C5 
very effectively. In system 1 the difference is more pronounced, as the non-bonded interactions 
between the methyl hydrogens and the methylene hydrogen are switched off (1 — 5 interaction, 
cf. Section 2). Comparison with experiment does not allow us to decide which scheme for the 
exclusion of nonbonded interactions is more appropriate. 

4.3 Segmental motion 

The double bonds reflect the reorientation of local segments. The ends reorient faster than inner 
monomers (Figure 12a); the six innermost monomers are comparable; chains of ten monomer 
length have a "bulk" inner part. 

There is a two step reorientation. On very short time scales there is a fast drop due to bond 
angle vibrations. This is not resolved here. For the inner monomers, a long-time process on the 
order of the reorientation of the whole chain (a few nano-seconds) sets in afterwards. 

There is little difference in the dynamics of chains of different lengths, at least in the limited 
range under study here: the end monomers are free to move regardless of the rest of the chain 
(Figure 12 b) and, for central monomers, the relaxation is the same within the (large) statistical 
error (Figure 12c). As seen in Figure 1 there are very few chains of any length other than 10 in 
system 2 and 3. 

Denault et al. estimated for chains of molecular weights between 7 x 10 3 g/mol and 1.5 x 
10 5 g/mol a segmental reorientation time of 1.0 ns at 303.15 K and of 43 ps at 373.15 K by 
analyzing their methylene reorientations. 9 For this, they used the Schaefer model 38 for segmental 
reorientation, where a ^-distribution Q f relaxation times is assumed, arising from cooperative 
local motion. The values are of the order of magnitude of the reorientation data presented here 
for 300 K. However, the chains are clearly longer and the focus is on the cis conformer. The 
chain length dependence is weak. An Arrhenius plot of the segmental reorientation time of 
Denault et. al. in comparison to the simulated reorientation of the double bond shows a similar 
temperature dependence (Figure 13). The activation energies deduced are E ( f m) w 33 kJ/mol 

and E^ xpS) = 65 kJ/mol at low temperature and E^*^ = 19 kJ/mol at high temperature. The 
effective experimental activation energy taking only the lowest and extrapolated highest point 
into account arrives at £^ xp " ) 36 kJ/mol rather close to the simulation value. To decide if 
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Figure 12: Reorientation of double bonds in system 2 at 300 K: a) Only the chains with ten monomers; 
results shown for different monomer distances from chain end. b) Only the first monomer for different 
chain lengths, c) Central monomers of the chains of length ten monomers and more (the first number 
denotes the chain length, the second the monomer number) in semilogarithmic representation. 
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Figure 13: Arrhenius-plot comparing the C=C reorientation time of the simulations (filled circles) 
to the segmental reorientation time inferred by Denault et al. from their experiments on cis-Pl at 
temperatures between 293 K and 373 K at molecular weights of 7000 to 130000. 9 For the simulations 
the values of table 8 are averaged at 300 K and 413 K respectively. The solid lines are exponential 
fits to the curves. 

there is a similar behavior as in experiments with two temperature regimes more simulations at 
intermediate temperatures would be necessary. 

5 Mapping onto a simple model 

As atomistic models such as the one presented here are obviously too demanding in computer 
time for long chain and/or long time investigations, mapping onto simpler models is highly 
desirable. The simple model we choose here is made of soft spheres (repulsive Lennard-Jones 
potential) connected by anharmonic springs and a weak harmonic bond angle potential leading 
to a persistence length of 1.5 monomeric units chosen similar to the atomistic model. For details 
of the model and results on structure and dynamics, see refs. 30,31. 

We have to choose a mapping of one model onto the other before we can compare. For the 
mapping of length scales the natural choice is the mean end-to-end distance. Thus, melts of 
10-mers of the atomistic and the simple models are simulated, then the lengths are set equal. 
The atomistic simulations of system 3b and simulations presented in Ref. 31 are used. The 
convergence of global chain properties in the atomistic case is shown in Figures 9b and 14. For 
the mapping of time scales a dynamic quantity is necessary. We map in Figure 14 the center- 
of-mass mean-square displacements g$ (t) onto each other. Due to limited simulation times, a 
corresponding mapping at 300 K could not be accomplished. With this mapping fixed, compar- 
ative analyses of the two models can be employed in order to check whether both models follow 
the same dynamics. Figure 15 shows the respective results for the mean-square-displacement of 
inner monomers, gi(t), the first three Rouse modes and the reorientation of nearest neighbor 
monomer-connecting vectors. These vectors are for the simple model just vectors connecting 
the beads. For the detailed model, these are vectors connecting the same carbon in adjacent 
monomers, i.e. Ci— Ci (dotted line) and C2— C2 (dashed line). For the Rouse mode analysis 
of the atomistic model the centers-of-mass of the double bond are taken as monomer positions. 
The first three Rouse modes and the inner-monomer-MSD do not coincide perfectly within this 
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Figure 14: Mean-square displacements of the center-of-mass in the atomistic simulation at 413 K 
(dashed line) and in the simple model with persistence length 1.5 monomer units (solid line). The 
ordinates are rescaled by the respective mean-square end-to-end distances. The time scale of the 
simple model is adjusted to bring the curves into coincidence for a whole order of magnitude. 

mapping. Still, they differ by only about a factor of two to three. The reorientation of the 
monomer- monomer vector even coincides without any refinement. The two differently defined 
vectors in the atomistic models are almost indistinguishable. Thus, on this scale the atomistic 
details are already rather unimportant; they may be incorporated into an intrinsic chain stiffness. 
These results show that the local dynamics are not exactly the same for the atomistic and the 
simple (coarse-grained) model. There are striking similarities, however, which are much stronger 
than one might expect as the models are completely different. This finding confirms the concept 
that the atomistic details on the scale of more than a monomer play only the role of shaping a 
persistence length. With the simple model we have carried out several investigations, 7 ' 30 ' 31 ' 39 
especially of local reorientation. 

One might argue, that discrepancies of a factor of three are no small effects. However, keeping 
in mind that the two models are drastically different one would expect on first sight a much 
stronger discrepancy. We performed the extreme step from the smallest possible model without 
quantum effects to a model where the complete identity of the polymer is only regarded by 
means of its stiffness. The successful mapping to the atomistic model is a further (retrospective) 
validation of those results. A more detailed investigation of the mapping will be shown in a 
separate publication. 40 



6 Conclusions 

This simulation study is concerned with decamers of irans-polyisoprene, which are short in 
comparison to experimental systems. Yet, the simulations are able to describe correctly local 
structural and dynamical features of the polymer. This can be seen by comparing results like 
chain extension, structure functions and correlation times for C— H vectors obtained by NMR. 
Taking into account that the simulated and experimental systems are not always identical in 
composition, temperature and so on, the agreement is quite good. We have also recently studied 
the free volume properties of our melts and compared them to positronium annihilation data, 41 
and the agreement is again very good. We thus conclude that our all- atom model provides a 
faithful description of this polymer. 
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Figure 15: Comparative analysis of the simple model with persistence length l p = 1.5 and the atomistic 
simulation at T=413 K. a) Inner monomer mean-square displacement (dashed line: atomistic model, 
center-of-mass of double bond of central monomer, solid line: simple model), b) Rouse modes (solid 
lines: first mode, dashed: second mode, dot-dashed: third mode), and c) reorientation of nearest 
neighbor monomer connecting vectors of the atomistic model at 413 K compared to the simple model. 



With the model, we analyzed the packing of chains. The detailed analysis of several inter- 
chain radial distribution functions shows that monomers approach each other most strongly with 
the exposed methyl groups followed by the methylene groups, whereas the vinyl carbons tend to 
be less accessible. All specificity in the interaction is however limited to the first solution shell. 
This is also seen in the mutual orientation of tangent vectors of neighboring chains. Chains 
have the tendency to be parallel at the first neighbor distance. The orientational correlation 
between second neighbors is already very small. From comparing different choices for the tangent 
vectors it is also evident that the orientations of bond vectors (i.e. short vectors) show some 
structure arising from atomistic interactions, whereas the orientation between inter- monomer 
(long) vectors is less structured and already close to what is found for a generic bead-spring 
model. 30 ' 39 

The reorientation dynamics of bond (C— C and C— H) vectors typically follows a two-stage 
process, a fast (picosecond) relaxation due to local vibrations followed by a long-time reorienta- 
tion characteristic for the reorientation of the parent polymer segment. The relative contribution 
of both processes to the overall reorientation is in good agreement with estimates from NMR 
measurements. 

In order to study the non-local structure and dynamics of irans-p olyisoprene extensive simu- 
lations have been performed with a generic bead-spring model which was only augmented by an 
intrinsic bending stiffness. 7,31 We have successfully mapped this model onto the present atom- 
istic simulation. After matching length and time scales, all characteristic time constants then 
agree to within a factor of 2 — 3. This illustrates that it is possible to develop polymer models 
at different levels of detail and have a description of polymer dynamics which smoothly connects 
both time scales. Systematic protocols for mapping atomistic to coarse-grained models and back 
are therefore being developed in our laboratory also for other polymer systems. 42-45 

An interesting technical point is the comparison between the polymer samples equilibrated 
with end-bridging Monte Carlo (EBMC) and those without. The EBMC proves to be a useful 
tool in the initial equilibration of polymer melts as it offers a way to wander through phase space 
more efficiently. At a global structural level, the EBMC provides better equilibrated samples. 
However, it is quite surprising that strictly local properties like monomer packing or bond- vector 
reorientation do not seem to be affected appreciably by the way the sample is prepared. 
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